Testing derivaties for 1D MT problem.

Especially the rx.projectFieldsDeriv


In [1]:
import SimPEG as simpeg
import simpegEM as simpegem, simpegMT as simpegmt
from SimPEG.Utils import meshTensor
import numpy as np

In [2]:
simpegmt.FieldsMT.FieldsMT_1D


Out[2]:
simpegMT.FieldsMT.FieldsMT_1D

In [3]:
# Setup the problem
sigmaHalf = 1e-2
# Frequency
nFreq = 33
# freqs = np.logspace(3,-3,nFreq)
freqs = np.array([100])
# Make the mesh
ct = 5
air = meshTensor([(ct,25,1.3)])
# coreT0 = meshTensor([(ct,15,1.2)])
# coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
core = np.concatenate( (  np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
bot = meshTensor([(core[0],10,-1.3)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
# Change to use no air
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core))], x0=x0)
# Make the model
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[ m1d.gridCC > 0 ] = 1e-8

rxList = []
for rxType in ['z1dr','z1di']:
    rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))
# Source list
srcList =[]
tD = False
if tD:
    for freq in freqs:
        srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
else:
    for freq in freqs:
        srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
# Make the survey
survey = simpegmt.SurveyMT.SurveyMT(srcList)

# Set the problem
problem = simpegmt.ProblemMT1D.eForm_psField(m1d)
problem.sigmaPrimary = sigma
problem.pair(survey)

# Get the fields
fields = problem.fields(sigma)

# Project the data
data = survey.projectFields(fields)

In [4]:
%debug


ERROR: No traceback has been produced, nothing to debug.

We need calculate this derivative. \begin{align} \underbrace{\frac{\partial P(f(u(m)),m^{fix})}{\partial f}}_{Rx} \end{align}

Use the rule \begin{align} \frac{d}{dx}\left( \frac{a(x)}{b(x)} \right) = \frac{\frac{d }{dx} a(x) b(x) - a(x)\frac{d }{dx} b(x) }{ b(x)^2 } \end{align}

In the case of the 1D MT problem the data is calculated as \begin{align} MT1Ddata = P(f(m)) &= \frac{P_{ex} f_e(src,m)}{P_{bx} f_b(src,m) \frac{1}{\mu_0}} = \frac{P_e u}{P_b f_b(u)} \\ \frac{\partial P(f(m))}{\partial u} v &= \frac{P_e}{P_b \frac{1}{mu_0} f_b(u)}v - \frac{P_e u}{\left(P_b \frac{1}{mu_0} f_b(u)\right)^2} P_b \frac{1}{mu_0} \frac{d f_b}{du} v \end{align} where u is the fields that we solve for. \begin{align} \frac{d f_b}{du} = - \frac{1}{i \omega} \nabla \end{align}


In [5]:
# Unused code &= \frac{ P_{ex} P_{bx} \frac{1}{\mu_0} \left( f_b(src,m) - f_e(src,m) \right) } { \left(P_{bx}f_b(src,m) \frac{1}{\mu_0} \right)^2 }

As matrices the formulas above can be written as \begin{align} \left[ \frac{\partial P(f(m))}{\partial u} v \right] = \left[ diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right] [P_e v] , diag[P_e u] diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right]^T diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right] \left[ P_b \frac{1}{mu_0} \frac{d f_b}{du}(v) \right] \right] \end{align}

The adjoint problem is done simliarly \begin{align} \left[ \frac{\partial P(f(m))}{\partial u} v \right]^T = [P_e]^T diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right]^T v - \left[ P_b \frac{d f_b}{du} \frac{1}{mu_0} \right]^T diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right] diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right]^T diag \left[ P_e u \right]^T v \end{align}


In [6]:
# def projectFields(self, src, mesh, u):
#         '''
#         Project the fields and return the
#         '''

#         if self.projType is 'Z1D':
#             Pex = mesh.getInterpolationMat(self.locs,'Fx')
#             Pbx = mesh.getInterpolationMat(self.locs,'Ex')
#             ex = Pex*mkvc(u[src,'e_1d'],2)
#             bx = Pbx*mkvc(u[src,'b_1d'],2)/mu_0
#             f_part_complex = ex/bx
#         real_or_imag = self.projComp
#         f_part = getattr(f_part_complex, real_or_imag)
#         return f_part

In [7]:
# Initate things for the derivs Test
src = survey.srcList[0]
rx = src.rxList[0]
v = np.random.randn(m1d.nN)
v0 = np.random.randn(m1d.nF+m1d.nE)
u0 = np.random.randn(m1d.nN)+np.random.randn(m1d.nN)*1j
f0 = problem.fieldsPair(m1d,survey)
f0[src,'e_1dSolution'] = u0
# f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0

In [8]:
# Run a test
def fun(u):
    f = problem.fieldsPair(m1d,survey)
    f[src,'e_1dSolution'] = u
    return rx.projectFields(src,m1d,f), lambda t: rx.projectFieldsDeriv(src,m1d,f0,t)
simpeg.Tests.checkDerivative(fun,u0,num=5,plotIt=False)


==================== checkDerivative ====================
iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order
---------------------------------------------------------
 0   1.00e-01    3.089e-04     2.730e-05      nan
 1   1.00e-02    2.841e-05     2.582e-07      2.024
 2   1.00e-03    2.818e-06     2.568e-09      2.002
 3   1.00e-04    2.816e-07     2.567e-11      2.000
 4   1.00e-05    2.816e-08     2.567e-13      2.000
========================= PASS! =========================
The test be workin!

Out[8]:
True

In [9]:
rx.projectFieldsDeriv(src,m1d,f0,u0)


Out[9]:
array([-0.00285083])

In [10]:
rx.projectFields(src,m1d,f0)


Out[10]:
array([[-0.00130618]])

In [11]:
# Test the Jvec derivative.

In [12]:
# print '%s formulation - %s' % (fdemType, comp)
CONDUCTIVITY = 0.01
m0 = np.log(np.ones(problem.mesh.nC)*CONDUCTIVITY)
# mu = np.log(np.ones(problem.mesh.nC)*MU)

if True:
    m0  = m0 + np.random.randn(problem.mesh.nC)*CONDUCTIVITY*1e-1 
#     mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1

# prb.mu = mu
# survey = prb.survey
def fun(x):
    
    return survey.dpred(x), lambda x: problem.Jvec(m0, x)
simpeg.Tests.checkDerivative(fun, m0, num=4, plotIt=False)


==================== checkDerivative ====================
iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order
---------------------------------------------------------
 0   1.00e-01    1.537e-05     1.005e-07      nan
 1   1.00e-02    1.541e-06     1.008e-09      1.999
 2   1.00e-03    1.541e-07     1.008e-11      2.000
 3   1.00e-04    1.541e-08     1.008e-13      2.000
========================= PASS! =========================
The test be workin!

Out[12]:
True

In [13]:
%debug


ERROR: No traceback has been produced, nothing to debug.

In [14]:
### Adjoint test

We have \begin{align} Jvec =&


In [15]:
# Run a test
TOL = 1e-4
FLR = 1e-20

def projectFieldsAdjointTest(fdemType, comp):
    print 'Adjoint %s formulation - %s' % (fdemType, comp)

    m  = np.log(np.ones(problem.mesh.nC)*0.01)
    if True:
        m  = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 

    u = problem.fields(m)
    v = np.random.randn(1)#+np.random.randn(1)*1j
    # print prb.PropMap.PropModel.nP
    w = np.random.randn(m1d.nN)+np.random.randn(m1d.nN)*1j

    vJw = v.dot(rx.projectFieldsDeriv(src,m1d,f0,w))
    wJtv = w.dot(rx.projectFieldsDeriv(src,m1d,f0,v,adjoint=True)).real
    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) 
    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
    return np.abs(vJw - wJtv) < tol
projectFieldsAdjointTest('e','projectFieldsDeriv')


Adjoint e formulation - projectFieldsDeriv
0.000177264079929 0.000177264079929 3.52365706058e-19 1e-07 True
Out[15]:
True

In [21]:
# Run a test
TOL = 1e-4
FLR = 1e-20

def getADeriv_mAdjointTest():
    print 'Adjoint test e formulation - getADeriv_m' 

    m  = np.log(np.ones(problem.mesh.nC)*0.01)
    if True:
        m  = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 

    u = problem.fields(m)
    v = np.random.randn(m1d.nN)#+np.random.randn(1)*1j
    # print prb.PropMap.PropModel.nP
    w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j

    vJw = v.dot(problem.getADeriv_m(freq,f0,w))
    wJtv = w.dot(problem.getADeriv_m(freq,f0,v,adjoint=True))
    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) 
    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
    return np.abs(vJw - wJtv) < tol
getADeriv_mAdjointTest()


Adjoint test e formulation - getADeriv_m
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-21-4c7e7058b8e9> in <module>()
     20     print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
     21     return np.abs(vJw - wJtv) < tol
---> 22 getADeriv_mAdjointTest()

<ipython-input-21-4c7e7058b8e9> in getADeriv_mAdjointTest()
     15     w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j
     16 
---> 17     vJw = v.dot(problem.getADeriv_m(freq,f0,w))
     18     wJtv = w.dot(problem.getADeriv_m(freq,f0,v,adjoint=True))
     19     tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])

/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.py in getADeriv_m(self, freq, u, v, adjoint)
     67         MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
     68         #
---> 69         u_src = u['e_1dSolution']
     70         dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
     71         if adjoint:

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Fields.pyc in __getitem__(self, key)
    150 
    151     def __getitem__(self, key):
--> 152         ind, name = self._indexAndNameFromKey(key, 'get')
    153         if name is None:
    154             out = {}

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Fields.pyc in _indexAndNameFromKey(self, key, accessType)
    131         srcTestList, name = key
    132         name = self._nameIndex(name, accessType)
--> 133         ind = self._srcIndex(srcTestList)
    134         return ind, name
    135 

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Fields.pyc in _srcIndex(self, srcTestList)
    100             ind = srcTestList
    101         else:
--> 102             ind = self.survey.getSourceIndex(srcTestList)
    103         return ind
    104 

/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Survey.pyc in getSourceIndex(self, sources)
    339         for src in sources:
    340             if getattr(src,'uid',None) is None:
--> 341                 raise KeyError('Source does not have a uid: %s'%str(src))
    342         inds = map(lambda src: self._sourceOrder.get(src.uid, None), sources)
    343         if None in inds:

KeyError: 'Source does not have a uid: e_1dSolution'

In [17]:
%debug


> /media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.py(69)getADeriv_m()
     68         #
---> 69         u_src = u['e_1dSolution']
     70         dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv

ipdb> u
> <ipython-input-16-0d25a7e16094>(17)getADeriv_mAdjointTest()
     16 
---> 17     vJw = v.dot(problem.getADeriv_m(freq,u0,w))
     18     wJtv = w.dot(problem.getADeriv_m(freq,u0,v,adjoint=True))

ipdb> d
> /media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.py(69)getADeriv_m()
     68         #
---> 69         u_src = u['e_1dSolution']
     70         dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv

ipdb> p u
array([  7.09826705e-01 -3.42993400e-01j,
         2.06495724e-02 -1.13236591e+00j,
         7.19940065e-01 +6.18299116e-01j,
        -3.76573540e-01 +9.09687373e-01j,
        -5.96166731e-01 +2.23309643e-01j,
         3.30533731e-01 +7.98068064e-01j,
        -1.41356795e+00 -1.74137927e-01j,
         1.51554175e+00 -2.16363375e-01j,
         1.59469627e+00 -1.87971566e+00j,
        -5.54637698e-01 -1.01936458e+00j,
         1.01869865e+00 +8.33443055e-01j,
        -4.92616671e-01 -1.97036224e+00j,
         1.00389102e+00 -1.52567261e-01j,
         7.90254779e-01 -5.68065046e-01j,
         9.17329964e-01 -1.18879097e+00j,
        -1.38317362e+00 -2.17789273e+00j,
         2.58324368e-01 -1.51003598e+00j,
        -2.04803634e+00 -1.44299943e-01j,
         3.33011285e+00 +1.52905384e+00j,
         2.91429972e-01 -1.91849967e+00j,
        -1.10273784e+00 -8.65794290e-01j,
        -2.80103156e-01 -1.42653192e+00j,
        -9.81404885e-01 +9.06845144e-02j,
         6.37579045e-01 +3.84148006e-01j,
        -6.83597988e-01 +6.91187965e-01j,
        -8.98575981e-01 +3.80347476e-02j,
         1.02647319e+00 -1.13876879e-01j,
         8.28006255e-01 -9.95903406e-01j,
         3.53931083e-01 -2.80188697e-01j,
        -1.09339690e+00 +1.43646692e+00j,
         2.59276626e+00 +4.09662561e-01j,
        -5.97873868e-01 -5.23852065e-01j,
         8.87054287e-01 -1.69872744e+00j,
         2.29767481e-01 -5.88143928e-01j,
        -6.97610315e-01 +5.68810360e-01j,
        -4.47049016e-01 +3.29354463e-01j,
        -3.31498311e-01 -8.27120276e-01j,
         1.57689998e-01 +6.68378872e-01j,
         9.47018631e-02 +2.60180446e-01j,
        -4.59177257e-01 -7.66896049e-01j,
        -1.16791204e+00 +4.20585804e-01j,
         2.72405147e-01 +3.48079787e-01j,
        -2.45880514e-02 +1.90226778e-01j,
        -7.03428363e-01 +1.20531568e+00j,
         1.60362217e+00 +6.88472259e-01j,
        -5.49144335e-01 +6.96208292e-01j,
        -9.83200346e-01 -1.15726575e+00j,
         1.22826524e+00 -8.58070855e-01j,
        -5.99020439e-01 -2.19182124e+00j,
        -5.69992871e-01 +1.68586900e+00j,
         9.26752644e-01 +1.05767829e+00j,
         8.64247444e-01 +1.51282356e-01j,
        -6.24547862e-01 -1.54517937e+00j,
        -1.72043526e+00 +1.91659834e+00j,
        -1.84501951e+00 -2.40165623e-01j,
         1.15210121e+00 -1.19627222e+00j,
        -3.93036395e-01 -1.94918396e+00j,
        -1.18266273e+00 -5.24652722e-01j,
         2.67461628e-01 -8.07279138e-01j,
        -1.56363641e+00 +1.59257615e-03j,
        -3.17211883e-01 -3.23379690e+00j,
        -1.08658939e-01 -8.66079039e-01j,
        -9.30350013e-01 -1.22668219e+00j,
        -1.07263137e+00 -5.56847000e-02j,
         1.27878013e+00 -1.34173537e+00j,
         7.91441492e-01 +6.42111414e-01j,
        -2.60121575e+00 -1.08140886e+00j,
        -1.15027808e+00 -4.47110142e-01j,
        -5.69747800e-01 +2.34065916e-01j,
         1.52155089e+00 -5.80340775e-01j,
         2.50052750e-01 -6.58793179e-01j,
         2.96944135e-01 +1.27455206e+00j,
         9.18939766e-01 -4.58122761e-01j,
         4.90015531e-01 -8.86288179e-01j,
        -1.65155083e+00 -6.57733967e-01j,
         9.72148169e-01 -1.19446820e+00j,
         2.08365350e+00 -1.53585595e-01j,
        -4.11101072e-01 +1.01418536e+00j,
         9.72580562e-01 -2.20470840e-01j,
         1.72922147e-01 -6.86360934e-01j,
         3.80255034e-01 +1.96862619e-01j,
        -1.53628572e+00 +5.31946011e-01j,
        -4.65939167e-01 +9.11725611e-01j,
        -1.06390365e+00 -3.07658816e-02j,
        -2.40734296e+00 -7.60530882e-01j,
         1.28242272e+00 +3.32085973e+00j,
         3.17004751e+00 +1.62654811e+00j,
        -8.49758476e-01 +9.86658359e-01j,
         4.99578557e-01 +4.01392488e-01j,
        -9.48695998e-01 +4.04322381e-01j,
         1.75371116e-01 +7.85458813e-01j,
         4.90232886e-01 -1.39280042e-01j,
        -2.31430252e-01 -8.44065839e-01j,
        -1.02319133e+00 +1.90560784e+00j,
         1.08940667e+00 +1.63046186e+00j,
        -4.34233324e-01 +2.90439312e-01j,
         4.47622371e-01 +1.26507976e+00j,
        -2.07317470e+00 -2.21558532e+00j,
        -3.19553487e-01 -9.19210427e-01j,
         6.09870390e-01 +3.05692898e-01j,
         4.92081860e-01 -2.18326059e-01j,
         1.48544627e+00 -2.02054490e-01j,
        -1.59612367e-01 +7.90853777e-02j,
         2.85854411e-01 +2.16567814e-01j,
         5.28063071e-01 +2.08838949e+00j,
         6.17580558e-01 +9.90748539e-01j,
        -2.64539211e-01 +1.16598665e+00j,
        -1.07357542e+00 -7.58002064e-01j,
         6.13399713e-02 +2.52338250e+00j,
        -4.84180135e-01 -1.45046020e+00j,
        -6.57519374e-01 -1.52719059e+00j,
        -1.54949084e-01 -7.64796563e-01j,
        -2.01172915e-02 +1.07377565e+00j,
         2.19840636e-01 +6.89622678e-01j,
        -2.64578500e+00 +1.46933966e+00j,
         1.22168696e+00 +2.33458454e-01j,
         2.71303177e+00 +2.77858685e-01j,
         4.39426665e-01 +1.79924695e-01j,
        -1.73858667e+00 -3.03473636e-01j,
         2.68736300e-01 +8.32359634e-01j,
         2.57517784e+00 +1.88097576e-01j,
         2.76775469e-01 -6.49891901e-01j,
        -9.87790301e-01 +4.50328269e-01j,
         9.19033669e-01 +1.01453748e+00j,
        -1.35706747e+00 +5.90240590e-01j,
         6.60118738e-01 +5.10361208e-01j,
        -2.17652642e-01 -4.26632365e-01j,
        -1.39163540e+00 +9.12739752e-01j,
        -1.76241829e+00 -1.04835908e+00j,
        -7.28279529e-01 +2.75961750e-01j,
        -1.32043513e+00 -3.28790183e-01j,
         5.88227093e-01 +9.41865970e-01j,
         1.97695607e+00 +1.04535450e+00j,
        -3.02599128e+00 +1.02031744e+00j,
        -1.34308763e+00 -4.48222977e-01j,
        -1.86687861e+00 +2.26767568e+00j,
         5.61472023e-01 -5.28652462e-01j,
        -5.26747946e-01 +2.75545054e-01j,
         6.16116514e-01 +3.41765793e-01j,
         1.38709750e+00 +5.87537003e-01j,
        -2.27923203e-01 +7.19317677e-01j,
        -8.66404713e-01 -1.80932081e-01j,
        -4.58146668e-01 +1.75955171e+00j,
        -5.12526647e-01 -1.91939928e-01j,
        -6.26532699e-01 -5.71301355e-01j,
        -3.35566040e-01 -5.06173544e-01j,
        -3.86124497e-01 -8.00560567e-01j,
         2.47252265e+00 +3.83928843e-01j,
         3.35561169e-01 +1.19309445e+00j,
        -6.19075444e-01 -1.19021628e+00j,
        -1.17335297e+00 +5.80500783e-01j,
        -2.35111248e-01 -1.59785495e+00j,
        -1.24819583e-01 -8.05281256e-01j,
        -9.18606131e-01 -1.25312332e-01j,
         3.14400972e-01 +1.18650079e+00j,
         6.50630049e-01 -6.99853497e-01j,
         1.28269722e+00 -1.29832850e+00j,
         2.07880202e-01 +2.56385692e-02j,
        -1.45089461e+00 -9.35093103e-01j,
         6.84146951e-01 +7.74422758e-01j,
         9.45439196e-01 -3.50695577e-01j,
         1.00057580e+00 -1.27412936e+00j,
        -2.01665427e-01 +1.73591867e+00j,
         1.76623768e+00 -1.42964638e+00j,
         1.12282181e+00 +9.02767375e-01j,
        -2.22663798e+00 +7.34672692e-02j,
         8.26230180e-02 -7.46601441e-01j,
        -4.92810641e-01 -5.49544995e-01j,
         1.06703148e-01 +1.70461156e+00j,
        -6.84949746e-01 -1.46950185e-01j,
         1.47509281e+00 +8.77421590e-01j,
        -9.83722735e-03 -6.30597166e-01j,
         2.36006569e-03 +2.12883413e-01j,
        -3.67899785e-01 -9.95508241e-02j,
         1.26752410e+00 +2.49370098e-01j,
        -1.19439605e+00 +2.94935402e-01j,
        -2.28617888e+00 +5.70564449e-02j,
         3.95617530e-02 -2.80995512e-01j,
        -7.51027805e-01 -8.14567233e-01j,
        -1.24160227e+00 -4.13553574e-01j,  -4.76626111e-01 +4.16480574e-01j])
ipdb> u_src
*** NameError: name 'u_src' is not defined
ipdb> c

In [18]:
# Run a test
TOL = 1e-4
FLR = 1e-20

def getRHSDeriv_mAdjointTest():
    print 'Adjoint test e formulation - getRHSDeriv_m'

    m  = np.log(np.ones(problem.mesh.nC)*0.01)
    if True:
        m  = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 

    u = problem.fields(m)
    v = np.random.randn(m1d.nN)#+np.random.randn(1)*1j
    # print prb.PropMap.PropModel.nP
    w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j

    vJw = v.dot(problem.getRHSDeriv_m(freq,w))
    wJtv = w.dot(problem.getRHSDeriv_m(freq,v,adjoint=True))
    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) 
    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
    return np.abs(vJw - wJtv) < tol
getRHSDeriv_mAdjointTest( )


Adjoint test e formulation - getRHSDeriv_m
(-10077.7224119-28916.4723594j) (-10077.7224119-28916.4723594j) (1.81898940355e-12-3.63797880709e-12j) 1.0 True
Out[18]:
True

In [ ]:
simpeg.mkvc(np.random.randn(survey.nD)+np.random.randn(survey.nD)*1j,2)

print survey.nD

In [19]:
TOL = 1e-4
FLR = 1e-20

def JvecAdjointTest():
    print 'Adjoint e formulation - Jvec' 

    m  = np.log(np.ones(problem.mesh.nC)*0.01)
    if True:
        m  = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 

    u = problem.fields(m)

    v = np.random.rand(survey.nD)
    # print prb.PropMap.PropModel.nP
    w = np.random.rand(problem.mesh.nC)

    vJw = v.dot(problem.Jvec(m, w, u))
    wJtv = w.dot(problem.Jtvec(m, v, u))
    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) 
    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
    return np.abs(vJw - wJtv) < tol

In [20]:
JvecAdjointTest()


Adjoint e formulation - Jvec
1.09508355274e-05 1.09508355274e-05 -1.01643953671e-20 1e-08 True
Out[20]:
True

In [ ]: